/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Declaration of the functions used to interface with Arepo.  Note
    that if you actually want to *use* the Arepo datastructures, you
    need to include the Arepo header files proto.h and allvars.h. We
    don't automatically do that here because they suck in a bunch of
    stuff like MPI that we don't want. */

// $Id$

#ifndef __mcrx_arepo__
#define __mcrx_arepo__


#include <utility>
#include "preferences.h"
#include <string>
#include "config.h"
#include "ray.h"
#include "mcrx-units.h"
#include "boost/thread/mutex.hpp"
#include <stdexcept>
#include <vector>

/** The arepo namespace contains functions that interface with the
    Arepo data structures. */
namespace arepo {
  void init(int*, char***);
  void load_snapshot(const std::string& snapfn, const std::string& parfile, 
		     int n_threads, const mcrx::T_unit_map& units);
  void calculate_midpoints();
  int locate (const mcrx::vec3d&, bool in_box=false);
  int locate_brute(const mcrx::vec3d& p);
  bool in_arepo_box(const mcrx::vec3d& p);
  int find_next_cell(const mcrx::ray_base& r, int p, int previous, 
		     mcrx::T_float& len);
  int find_closest_neighbor_in_cell(const mcrx::vec3d& p, int pp);
  mcrx::vec3d extend_to_box(const mcrx::ray_base& r);
  mcrx::vec3d get_point(int p);
  bool primary_cell(int dp_idx);
  int sph_index(int dp_idx);
  mcrx::T_float column_to_distance_fraction(int dp_idx, mcrx::T_float fn, 
					    const mcrx::vec3d& pos, 
					    const mcrx::vec3d& dir,
					    mcrx::T_float L);
  mcrx::T_float cell_density(int sphp_idx, const mcrx::vec3d& p);

#ifndef NDEBUG
  void assert_contains(int dp_idx, mcrx::vec3d p);
#else
  // Assert_contains is inlined to be a nop if ndebug is set.
  inline void assert_contains(int dp_idx, mcrx::vec3d p) {};
#endif

  class arepo_pressure_functor;
  void ensure_positive_densities();

  // globals for unit conversion
  extern mcrx::T_unit_map arepo_units;
  extern mcrx::T_float lcon;
  extern mcrx::T_float mcon;
  extern mcrx::T_float tcon;

  // Voronoi cell face data structures
  extern std::vector<int> primary_cells;
  extern std::vector<std::pair<int, int> > midpoint_idx;
  extern std::vector<mcrx::vec3d> midpoints;
  extern std::vector<int> opposite_points;

  // Arepo box extent in Sunrise units
  extern mcrx::vec3d arepomin, arepomax;
};


class arepo::arepo_pressure_functor {
  const mcrx::T_float pcon, dcon;

public:
  arepo_pressure_functor();
  arepo_pressure_functor(const arepo_pressure_functor& rhs) :
    pcon(rhs.pcon), dcon(rhs.dcon) {};
  ~arepo_pressure_functor() {};

  std::pair<mcrx::T_float,std::pair<mcrx::T_float, mcrx::T_float> > 
  operator()(const mcrx::vec3d& pos) const;
};

#endif
